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Abstract 


A  previously  developed  intermolecular  potential  for  nitramines  and  several  other  classes  of 
nitrocompound  crystals  has  been  used  to  investigate  the  behavior  of  the  energetic  materials 
hexahydro- 1,3,5  -trinitro- 1,3,5  -5-trazine  (RDX)  ,1,3,5 ,7-tetranitro- 1 ,3 ,5 ,7 -tetraazacyclo-octane 
(HMX),  2,4,6,8,10,12-hexanitrohexaazaisowurtzitane  (HNIW),  and  pentaerythritol  tetranitrate 
(PETN)  under  hydrostatic  compression.  Isothermal-isobaric  molecular  simulations  (assuming 
the  rigid-molecule  approximation)  molecular-packing  calculations  were  used  to  perform  the 
analyses.  In  the  case  of  the  RDX,  HMX,  and  HNIW  crystals,  the  results  indicate  that  the 
proposed  potential  model  is  able  to  accurately  reproduce  the  changes  in  the  structural 
crystallographic  parameters  as  functions  of  pressure  for  the  entire  range  of  pressures  that  has 
been  investigated  experimentally.  In  addition,  the  calculated  bulk  moduli  of  RDX  and  HMX 
were  found  to  be  in  good  agreement  with  the  corresponding  experimental  results.  In  the  case  of 
the  PETN  crystal,  the  crystallographic  parameters  have  been  reproduced  with  an  acceptable 
accuracy  at  pressures  up  to  about  5  GPa.  The  larger  deviations  from  the  experimental  results  at 
greater  pressures  indicate  the  limitations  of  the  rigid-molecule  model  when  applied  to  floppy 
molecules.  The  similarity  of  the  results  determined  in  molecular-packing  calculations  relative 
to  those  from  molecular  dynamics  simulations  suggest  that  the  former  method  can  be  used  as  an 
efficient  tool  for  rapid  tests  of  the  crystal  structure  modification  under  pressure. 
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1.  Introduction 


New  strategies  for  the  development  of  new  energetic  materials  or  deployment  of  existing 
materials  in  advanced  weapons  platforms  have  incorporated  technologies  that  result  in  cost  and 
time  efficiency.  These  strategies  include  modeling  and  simulation  at  the  various  stages  in  the 
developmental  process.  Modeling  and  simulation  are  recognized  to  be  a  cost-effective  and 
efficient  means  of  achieving  such  goals  in  any  developmental  process.  Atomistic  simulation  for 
characterization  and  prediction  of  physical  and  chemical  behavior  of  energetic  materials 
promises  to  be  one  of  the  more  powerful  and  effective  modeling  methodologies  that  is 
incorporated  into  the  new  developmental  strategies.  Accurate  atomistic  predictions  provide 
information  of  the  fundamental  mechanisms  of  processes  that  determine  performance  of  the 
materials.  However,  the  effectiveness  of  the  simulations  is  limited  by  the  accuracy  of  the 
description  of  the  interaction  potential  of  the  models.  An  attempt  has  been  made  to  develop 
classical  models  of  energetic  materials  that  will  accurately  reproduce  known  properties  of  these 
materials.  The  approach  has  been  to  first  develop  simple  potential  functions  to  describe  the 
interactions  between  molecules  in  the  crystals,  while  the  overall  goal  is  to  enhance  these  models 
such  that  the  accurate  description  of  different  complex  physical  and  chemical  processes, 
including  chemical  reactions  in  the  condensed  phase,  can  be  achieved.  Since  most  of  the 
processes  of  interest  here  are  in  the  condensed  phase  and  involve  systems  containing  large 
polyatomic  molecules,  it  is  imperative  that  the  interactions  are  described  by  simple  functions; 
otherwise  atomistic  simulation  could  become  computationally  intractable.  Such  simple  functions 
have  been  developed  and  have  been  evaluated  in  a  series  of  studies  that  predict  crystallographic 
parameters  at  ambient  conditions  [1-5]. 

An  initial  study  [1]  has  shown  how  an  atom-atom  (6-exp)  Buckingham  potential  with 
electrostatic  interaction  terms  in  the  form  of  partial  charges  associated  with  the  atoms  of  the 
molecules  can  be  parameterized  to  reproduce  the  experimental  crystal  structure  of  the  a-form  of 
the  solid  explosive,  hexahydro-l,3,5,-trinitro-l,3,5-s-triazine  (RDX).  This  intermolecular 
potential  was  used  to  simulate  the  RDX  crystal  structure  in  isothermal-isobaric  molecular 
dynamics  (NPT-MD)  calculations  at  ambient  pressure  and  for  temperatures  ranging  from  4.2  to 
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325  K.  The  results  of  the  simulations  indicated  very  good  agreement  with  experiment,  with  the 
lattice  dimensions  being  within  2%  of  experiment  and  almost  no  rotational  or  translational 
disorder  of  the  molecules  in  the  unit  cell.  The  space-group  symmetry  was  maintained  throughout 
the  simulations  for  the  average  structures.  Additionally,  the  predicted  thermal  expansion 
coefficients  were  in  reasonable  agreement  with  experiment. 

The  utility  of  the  proposed  potential  was  expanded  when  it  was  shown  that  the  same 
Buckingham  6-exp  potential  terms  can  be  used  without  modification  to  characterize  (through 
molecular  packing  [MP]  and  NPT-MD  simulations)  the  structures  and  their  thermal  dependence 
for  other  nitramine  crystals  (i.e.,  2,4,6,8,10,12-hexanitrohexaazaisowurtzitane  [HNIW]  [2]  and 
l,3,5,7-tetranitro-l,3,5,7-tetraazacyclo-octane  [HMX]  [3]).  Investigations  indicate  that  this 
potential  predicts  accurately  not  only  the  crystallographic  structures  of  different  phases  of  these 
crystals  but  also  the  correct  lattice  energies  and  the  relative  order  of  stability.  Particularly,  the 
potential  indicate  the  stability  order  e  >  |3  >  y  and  P  >  a  >  5  for  the  corresponding  phases  of 
HNIW  and  HMX,  in  agreement  with  the  experimental  measurements  [6,  7]. 

More  recently,  investigations  of  the  transferability  of  this  interaction  potential  in  molecular 
simulations  of  30  nitramine  crystals  [4]  have  been  extended.  The  molecules  associated  with  the 
nitramine  crystals  were  chosen  as  representative  examples  of  acyclic  and  cyclic  nitramines.  In 
the  latter  case,  different  types  of  mono-  and  polycyclic  nitramines  have  been  included, 
particularly  crystals  of  importance  in  energetic  materials.  For  most  of  the  crystals,  the  predicted 
structural  lattice  parameters  differ  by  less  than  2%  from  the  experimental  structures,  with  small 
rotations  and  practically  no  translations  of  the  molecules  in  the  asymmetric  unit  cell. 

Finally,  studies  have  expanded  to  assess  whether  the  interaction  potential  could  be  used  to 
model  crystals  beyond  the  class  of  nitramines.  MP  calculations  have  been  performed  on 
51  crystals  comprising  a  wide  variety  of  compounds  such  as  nitroalkanes,  nitroaromatics, 
nitrocubanes,  polynitroadamantanes,  polynitropolycycloundecanes,  polynitropolycyclodode- 
canes,  hydroxy-nitroderivatives,  nitrobenzonitriles,  nitrobenzotriazoles,  and  nitrate  esters,  such 
that  a  comprehensive  test  to  this  potential  was  achieved  [5],  It  was  shown  that,  for  the  majority 
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of  these  crystals,  the  potential  model  accurately  reproduces  the  crystallographic  structural  and 
energetic  data  determined  experimentally.  Moreover,  in  the  same  study,  the  temperature 
dependence  of  the  crystallographic  parameters  has  been  analyzed  for  two  important  energetic 
crystals,  2,4,6-trinitrotoluene  (TNT)  in  the  monoclinic  phase  and  the  pentaerythritol  tetranitrate 
(PETN)  crystal  in  the  tetragonal  phase.  In  each  case,  the  results  show  that,  throughout  the  MD 
simulations,  the  average  structures  of  the  crystals  maintain  the  same  space  group  symmetry  as 
the  one  determined  experimentally  and  there  is  a  good  agreement  between  the  calculated 
crystallographic  parameters  and  the  experimental  values. 

The  present  paper  considers  another  category  of  tests  through  which  one  can  assess  the 
quality  of  the  potentials.  In  particular,  focus  is  on  the  analysis  of  the  hydrostatic  compression  of 
some  energetic  materials  and  if  the  structural  changes  observed  experimentally  can  be  described 
accurately  with  the  present  set  of  potentials.  For  this  purpose,  consideration  is  given  to  the  case 
of  the  nitramine  crystals  RDX  (a-phase),  HMX  (3-phase),  and  HNIW  (e-phase)  and  the 
non-nitramine  crystal  PETN  for  which  experimental  information  is  available.  The  configurations 
of  the  molecules  corresponding  to  these  are  illustrated  in  Figure  1.  As  can  be  seen,  RDX  and 
HMX  have  monocyclic  molecular  configurations,  HNIW  is  polycyclic,  and  PETN  is  acyclic. 

As  in  the  preceding  studies  [1-5],  the  previously  presented  [1]  RDX  Buckingham  potential 
was  used,  plus  Coulombic  interactions  terms  obtained  through  fitting  of  partial  charges  centered 
on  each  atom  of  the  molecule  (in  the  experimental  arrangement)  to  a  quantum  mechanically 
derived  electrostatic  potential  [8].  It  has  been  shown  for  both  the  nitramine  [4]  and 
non-nitramine  [5]  crystals  that  the  best  agreement  between  the  calculated  and  experimental 
energies  is  obtained  when  the  set  of  charges  is  determined  using  methods  that  employ  electron 
correlation  effects.  For  example,  in  the  case  of  a  set  of  30  nitramine  crystals  previously  analyzed 
[4],  an  average  deviation  of  the  Hartree-Fock  lattice  energies  of  12.8%  was  found  from  the 
corresponding  Moller-Plesset  (MP2)  [9-12]  energies.  The  use  of  density  functional  theory 
(B3LYP)  to  evaluate  the  electrostatic  charges  decreases  these  deviations  of  the  lattice  energies  to 
about  2.6%.  Thus,  in  the  present  case,  the  sets  of  charges  were  derived  using  the  second-order 
Moller-Plesset  (MP2)  perturbation  theory. 
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Figure  1.  Molecular  Configurations  of  (a)  RDX  (a-Phase),  (b)  HMX  (P-Phase),  (c)  HNIW 
(e-Phase),  and  (d)  PETN  (Tetragonal  Phase). 


As  in  the  previous  studies  [1-5],  the  main  limitation  of  these  models  remains  the  assumption 
of  rigid  molecules.  However,  as  previously  discussed  by  Pastine  [13],  the  initial  compression  of 
an  organic  explosive  is  almost  entirely  due  to  a  reduction  of  intermolecular  distances.  It  is  only 
at  small  intermolecular  distances  that  the  increase  of  van  der  Waals  repulsions  become 
comparable  to  the  intramolecular  repulsions  along  the  covalent  bonds  [13].  Consequently,  the 
assumption  of  rigid  molecules  should  be  adequate  in  describing  the  initial  compression  of  such 
crystals  within  the  regime  in  which  intramolecular  deformations  are  small.  This  should  be  valid 
for  the  crystals  containing  molecules  that  are  not  very  floppy,  such  as  the  RDX,  HMX,  and 
HNIW  crystals.  This  approach,  in  which  rigid  molecules  are  assumed  in  simulations  of  the 
hydrostatic  effects  on  crystallographic  parameters,  has  also  been  used  previously  [14]. 
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The  organization  of  the  paper  is  as  follows.  In  section  2,  the  intermolecular  potential  used  to 
simulate  the  crystals  is  presented.  In  section  3,  the  details  of  calculations  using  molecular 
packing  methods  and  isothermal-isobaric  MD  calculations  are  described.  The  results  of  these 
calculations  are  given  in  section  4.  The  main  conclusions  are  summarized  in  section  5. 

2.  Intermolecular  Potential 

The  intermolecular  potential  used  in  the  present  study  is  that  previously  used  for  nitramine 
[4]  and  non-nitramine  [5]  crystals.  Therefore,  only  brief  details  are  provided.  The 
intermolecular  interactions  between  the  molecules  of  the  crystal  are  approximated  by  using 
superpositions  of  pairwise  Buckingham  (6-exp)  (repulsion  and  dispersion)  and  Coulombic 
(C)  potentials  of  the  form 


V*(r)  = 


(1) 


and 


Q  a  p 
47t7tor  ’ 


(2) 


where  r  is  the  interatomic  distance  between  atoms  a  and  P,  q«  and  qp  are  the  electrostatic  charges 
on  the  atoms,  and  £o  is  the  dielectric  permittivity  constant  of  vacuum. 

The  parameters  for  the  6-exp  potential  in  equation  (1)  are  those  previously  determined  for  the 
RDX  crystal  [1].  The  heteroatom  parameters  are  calculated  from  the  homoatom  parameters 
using  the  same  combination  rules  as  previously  reported  [1]. 

The  assignments  of  the  electrostatic  charges  were  made  by  using  the  set  of  atom-centered 
monopole  charges  for  the  isolated  molecule  (with  the  structure  fixed  at  the  experimental 
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crystallographic  configuration)  that  best  reproduces  the  quantum  mechanically  derived 
electrostatic  potential.  The  electrostatic  potential  is  calculated  over  grid  points  surrounding  the 
van  der  Waals  surface  of  the  molecules.  This  method  of  fitting  the  electrostatic  potential  was 
proposed  by  Breneman  and  Wiberg  [8]  and  is  incorporated  in  the  Gaussian  94  package  of 
programs  [15]  under  the  keyword  CHELPG  (electrostatic-potential-derived  atomic  charges). 
These  calculations  have  been  done  at  MP2  theoretical  level  using  a  reasonable  quality  basis  set 
(i.e.,  6-31G**  [split-valence  plus  d-type  and  p-type  polarization  functions])  [16]. 

3.  Computational  Approach 

3.1  Molecular  Packing  Calculations.  Preliminary  tests  of  the  capability  of  the  interaction 
potential  to  adequately  predict  the  crystal  structural  behavior  under  hydrostatic  compression  are 
performed  using  molecular  packing  calculations  [17,  18],  in  which  the  lattice  energy  of  a  crystal 
is  minimized  with  respect  to  its  structural  degrees  of  freedom.  These  calculations  have  been 
done  using  the  algorithm  proposed  by  Gibson  and  Scheraga  [19]  for  efficient  minimization  of  the 
energy  of  a  fully  variable  lattice  composed  of  rigid  molecules  and  implemented  in  the  program 
LMIN  [20].  The  nonbonded  interactions  were  attenuated  using  a  cubic  spline  function  from  Pa 
to  Qa,  to  ensure  the  continuity  of  the  function  and  its  first  derivative.  Here  a  is  the  value  of  r  in 
equation  (1)  at  which  Vap(r)  =  0  and  dV«p(r)/dr<0.  The  parameters  P  and  Q,  which  specify  the 
start  and  the  end  of  the  cubic  feather  (see  Sorescu  et  al.  [1]  and  Desiraju  [18]  for  details)  were  set 
to  20.0  and  20.5,  respectively.  The  Coulombic  potential  terms  of  the  form  given  in  equation  (2) 
are  summed  over  the  lattice  using  the  Ewald  technique  as  previously  described  [1].  Finally,  the 
effect  of  pressure  on  the  crystallographic  parameters  has  been  simulated  by  adding  a  potential 
term  of  the  form  P(V  -  Vo)  [21],  where  Vo  is  the  volume  of  a  suitably  chosen  unit  cell  at  zero 
pressure. 


3.2  Constant  Pressure  and  Temperature  Molecular  Dynamics  Calculations.  A  more 
comprehensive  test  of  the  ability  of  the  interaction  potential  to  predict  the  crystal  structure  of  the 
molecular  crystals  under  hydrostatic  compression  has  been  done  using  constant  pressure  and 
temperature  (NPT)  molecular  dynamics  simulations,  in  which  there  are  no  geometric  constraints 
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other  than  the  assumption  of  rigid-body  molecules.  This  method  yields  average  equilibrium 
properties  of  the  lattice  as  functions  of  temperature  and  pressure. 

The  Nose-Hoover  barostat  algorithm  [22]  has  been  used  as  implemented  in  the  program 
DL_POLY_2.0,’  to  simulate  the  crystals  at  various  temperatures  and  pressures.  In  this  case,  the 
equations  of  motion  for  both  the  translation  of  rigid  molecules  and  the  simulation  cell  are 
integrated  using  the  Verlet  leap-frog  scheme  [23].  The  molecular  rotational  motion  is  handled 
using  Fincham’s  implicit  quaternion  algorithm  [24]. 

The  MD  simulation  cells  consist  of  boxes  containing  3x3x3,  4x2x3,  3x2x3,  and 
3x3x4  crystallographic  unit  cells  for  the  RDX,  HMX,  HNIW,  and  PETN  crystals, 
respectively.  These  choices  of  the  simulation  boxes  ensure  the  use  of  a  cutoff  distance  for  the 
intermolecular  potentials  of  about  10  A.  For  each  of  the  four  different  molecular  crystals,  a 
simulation  corresponding  to  the  lowest  pressure  was  first  performed,  with  the  position  and 
orientation  of  the  molecules  in  the  unit  cell  initially  set  to  be  identical  to  those  for  the 
experimental  structure.  The  systems  were  then  equilibrated  at  298  K  and  atmospheric  pressure. 
In  all  production  runs,  the  total  integration  time  corresponded  to  12,000  time  steps  (1  time 
step  =  2  x  10-15  s),  of  which  2,000  steps  were  equilibration.  The  velocities  were  scaled  after 
every  five  steps  during  the  equilibration  period  so  that  the  internal  temperature  of  the  crystal 
mimicked  the  imposed  external  temperature.  Properties  were  then  calculated  and  accumulated 
for  averaging  over  the  next  1,0000  integration  steps  in  the  simulation.  In  subsequent  runs, 
performed  at  successively  higher  pressures  and  a  constant  temperature  of  298  K,  the  initial 
configurations  of  the  molecular  positions  and  velocities  were  those  corresponding  to  the  final 
values  from  the  preceding  lower-pressure  simulation. 

The  lattice  sums  were  calculated  subject  to  the  use  of  minimum-image  periodic  boundary 
conditions  in  all  dimensions  [23].  The  interactions  were  determined  between  the  sites  (atoms)  in 
the  simulation  box  and  the  nearest-image  sites  within  the  cutoff  distance.  In  these  calculations, 
the  Coulombic  long-range  interaction  were  handled  using  Ewald’s  method  [23]. 

DLJPOLY  is  a  package  of  molecular  simulation  routines  written  by  W.  Smith  and  T.  R.  Forester,  copyright  The 
Council  for  the  Central  Laboratory  of  the  Research  Councils,  Daresbury  Laboratory  at  Daresbury,  Nr., 
Warrington,  1996. 
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The  main  quantities  obtained  from  these  simulations  were  the  average  lattice  dimensions  and 
the  corresponding  volume  of  the  unit  cell.  Additional  information  about  the  structure  of  the 
crystal  has  been  obtained  by  calculating  the  center-of-mass  (COM)-COM  radial  distribution 
functions  (RDF).  Such  quantities  have  been  calculated  from  recordings  done  at  every  10th  step 
during  trajectory  integrations. 

4.  Results  and  Discussions 

4.1  RDX  Crystal.  Crystalline  RDX  exists  in  two  phases  [25]:  the  ambient  phase  (a-solid), 
for  which  the  structure  has  been  characterized  by  neutron  diffraction  measurements  [26],  and  an 
unstable  phase  O-solid),  the  crystal  structure  of  which  has  not  been  determined.  The  structure  of 
ct-RDX  at  room  temperature  and  1-atm  pressure  belongs  to  the  orthorhombic  space  group  Pbca 
with  Z  =  8  molecules  per  unit  cell.  The  linear  and  volume  compression  of  RDX  have  been 
investigated  by  Olinger  et  al.  [27]  for  pressures  up  to  9  GPa  using  a  high-pressure  x-ray 
diffraction  technique.  It  has  shown  that  the  ambient  RDX  polymorph  is  stable  until  a  pressure  of 
3.95  GPa.  Above  this  pressure,  a  new  polymorph  phase  is  formed,  which  remains  stable  until 
9  GPa.  The  present  study  does  not  consider  this  change  to  the  new  polymorphic  state,  so  the 
range  of  pressures  investigated  is  up  to  3.95  GPa. 

The  results  of  MP  and  NPT-MD  calculations  are  summarized  in  Table  1  and  compared  with 
experimental  data  in  Figure  2.  At  P  =  0  the  relative  differences  between  the  predicted  and  the 
experimental  values  are  very  small.  By  considering  as  the  reference  for  comparison  the  results 
determined  by  Olinger  et  al.  [27],  the  percentage  errors  for  lattice  dimensions  a ,  b,  and  c  are 
0.64%,  0.47%,  and  -1.01%  for  MP  results  and  1.52%,  1.73%,  and  0.12%  for  the  MD-NPT  data. 
The  increase  of  pressure  from  0  to  3.95  GPa  does  not  significantly  change  the  differences 
between  the  predicted  and  the  experimental  sets  of  values.  For  example,  at  the  largest  pressure 
considered  here  (3.95  GPa),  the  corresponding  percent  deviations  are  0.66%,  2.23%,  and  0.29% 
for  the  MP  values  and  0.91%,  2.60%,  and  0.54%  for  the  MD-NPT  data.  It  appears  that  there  is  a 
slight  increase  of  the  deviation  from  the  experimental  data  with  pressure  for  the  b  lattice 
dimensions  only;  the  relative  differences  between  experiment  and  predictions  of  the  other  lattice 
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Table  1.  Lattice  Parameters  Obtained  in  Crystal  Packing  and  Molecular  Dynamics 
Calculations  for  the  a-RDX  Crystal  as  a  Function  of  Pressure  at  T  =  298  Ka 


Pressure 

(GPa) 

Molecular  Packing 

NPT-MD 

Lattice  Lengths 

Unit  Cell 
Volume 
(A3) 

Unit  Cell 
Volume 
(A3) 

a 

0 

(A) 

b 

o 

(A) 

c 

0 

(A) 

a 

0 

(A) 

b 

o 

(A) 

dUi 

13.182 

11.574 

10.709 

1633.856 

— 

— 

— 

— 

WiM 

13.202 

11.596 

10.717 

1640.670 

— 

— 

— 

— 

0.00 

13.2863 

11.6510 

10.6086 

11.7974 

1596.575 

1635.074 

pSillB 

13.1581 

11.4986 

10.4395 

1579.493 

1612.963 

1.75 

13.0005 

11.4283 

10.3277 

1540.740 

2.75 

12.8624 

11.2696 

10.1669 

1473.734 

1490.185 

3.36 

1451.978 

12.8399 

11.2638 

10.1361 

1465.962 

3.95 

12.7535 

11.1742 

10.0589 

11.2139 

10.0838 

1445.703 

aThe  experimental  values  indicated  correspond  to  T  =  298  K  and  atmospheric  pressure. 
b  Neutron-diffraction  values  from  Choi  and  Prince  [26]. 
c  X-ray  diffraction  values  from  Olinger  et  al.  [27]. 


dimensions  remain  practically  unchanged.  The  corresponding  parameters  of  the  quadratic  fits  of 
the  predicted  lattice  parameters  and  volume  as  function  of  pressure  are  given  in  Table  2.  An 
important  finding  is  that  lattice  dimensions  predicted  in  MP  calculations  are  very  close  to  the 
experimental  data.  This  is  notable  since  MP  calculations  require  significantly  less  computational 
time  than  the  corresponding  NPT-MD  calculations. 


The  normalized  dependence  of  the  unit  cell  volume  on  pressure  is  shown  in  Figure  2(b).  The 
volume  compressibility  predicted  in  from  NPT-MD  simulations  is  very  close  to  that  seen 
experimentally.  The  calculated  ratio  V/Vo  obtained  in  MD  simulations  at  3.95  GPa  is  0.852, 
while  the  corresponding  experimental  value  [27]  is  0.846.  At  this  pressure,  the  predicted  MP 
ratio  is  slightly  larger,  with  a  value  of  0.873. 


The  dependence  of  the  NPT-MD  unit  cell  volume  on  pressure  was  determined  by  fitting  the 
calculated  values  to  the  Mumaghan  equation  [28]: 
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Table  2.  Coefficients  of  the  Quadratic  Fits  of  the  Form  ao+aiP+a2P2  for  RDX,  HMX,  and 
HNIW  and  of  the  Cubic  Fit  of  the  Form  a0(l+aiP+a2P2+«3P3)  for  PETN  of  the 
Lattice  Constants  and  Unit  Cell  Volume  as  Function  of  Pressure  (GPa)  Using 
Results  From  the  NPT-MD  Calculations.3  The  Calculated  Bulk  Modulus  (B0),  its 
Pressure  Derivative  at  Zero  Pressure  (B0’)»  and  Zero  Volume  Coefficient  Using 
Equation  (3)  Are  Indicated  Together  With  the  Corresponding  Experimental 
Values  Where  Available 


System 

0o 

d\ 

a2 

Boexp 

(GPa) 

B'oexp 

RDX 

_ _ 1 

a 

13.3964 

-0.23490.0201 

0.0201 

— 

— 

— 

_ 

b 

11.7786 

-0.24620.0269 

0.0269 

— 

— 

— 

— 

1 

c 

10.7103 

-0.2699 

0.0289 

— 

— 

— 

— 

E9 

V 

1689.604 

7 

-104.9649 

11.1950 

1639.62 

12.93 

6.77 

13.0a 

P-HMX 
(P2]/c  setting) 

a 

6.5746 

-0.0947 

— 

— 

— 

b 

11.0206 

-0.1488 

■ 

— 

— 

— 

c 

9.0416 

-0.1252 

0.0078 

— 

— 

— 

V 

531.1252 

-22.72731 

1.5919 

534.88 

14.53 

9.57 

ESIB 

P-HMX 
(P2|/n  setting) 

a 

6.5069 

-0.0803 

0.0046 

— 

— 

— 

b 

10.9275 

-0.1317 

0.0070 

l 

— 

— 

_ 

c 

7.3745 

-0.1061 

0.0075 

— 

— 

— 

_ 

V 

516.3269 

-19.4548 

1.2220 

519.79 

16.86 

9.50 

— 

— 

e-HNIW  I 

a 

8.8956 

-0.1215 

0.0113 

— 

— 

— 

— 

_ 

b 

12.5774 

-0.2564 

0.0288 

— 

— 

— 

— 

_ 

c 

13.5481 

-0.2475 

0.0257 

— 

— 

— 

— 

_ 

V 

1461.6076 

-75.4732 

8.3433 

1463.99 

15.58 

9.37 

— 

— 

— 

— 

— 

1465.49 

14.67 

9.93 

— 

— 

PETN 

0o 

0i 

02 

03 

Vof,t 

(A3) 

Bo 

(GPa) 

B0' 

Boexp 

(GPa) 

B'oexp 

a  (md) 

9.3348 

1.552  xlO-2 

1.863  x  1  O'3 

-0.920  x  10'4 

— 

— 

_ 

_ 

a  (exp)b 

9.3830 

-2.052  x  10‘2 

2.230  x  10“3 

-1.041  x  10-4 

— 

— 

— 

— 

c  (md) 

6.6500 

-1.921  x  10‘2 

2.101  x  10"3 

-1.021  x  10-4 

— 

— 

— 

— 

c  (exp)b 

6.7150 

-2.832  x  10“2 

3.295  x  10'3 

-1.458  x  10’4 

— 

— 

— 

— 

V  (md) 

578.7628 

-4.8934  x  10~2 

5.839  x  10~3 

-2.880  x  10“4 

582.49 

14.09 

10.39 

9.9° 

11. oc 

a  The  calculated  bulk  modulus  (B0),  its  pressure  derivative  at  zero  pressure  (B0'),  and  zero  volume  coefficient  using 
equation  (3)  are  indicated  together  with  the  corresponding  experimental  values  where  available. 
bData  from  Olinger  et  al.  [27]. 
c  Data  from  Olinger  et  al.  [32]. 
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• 

a  (EXP) 

o 

a  (MD) 

T 

a  (MP) 

V 

b  (EXP) 

■ 

b  (MD) 

□ 

b  (MP) 

♦ 

c  (EXP) 

o 

c  (MD) 

▲ 

1 

c  (MP) 

Figure  2.  Comparison  of  the  Crystallographic  Parameters  Obtained  in  MP  and  NPT-MD 
Simulations  for  the  a-RDX  Crystal  With  the  Experimental  Results.  Dependence 
of  (a)  Lattice  Dimensions  and  (b)  Volume  Compression  V/Vo  on  the  External 
Pressure. 


(3) 


In  this  equation,  V  is  the  volume  at  pressure  P,  Vo  is  the  fitted  volume  at  P  =  0,  Bo  is  the  bulk 
modulus,  and  Bq  =dB0/dP.  The  best-fit  parameters  are  given  in  Table  2.  The  predicted  bulk 
modulus  and  its  pressure  derivative  are  very  close  to  the  corresponding  experimental  values  with 
relative  percentage  errors  of  -0.5%  and  2.5%,  respectively. 
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4.2  HMX  Crystal.  Crystalline  HMX  can  exist  in  four  polymorphic  phases,  known  as  the  a, 
P,  y,  and  8  forms  [3].  The  stable  form  at  room  temperature  is  (3-HMX.  It  has  a  monoclinic 
structure  with  P2,/n  symmetry  [29]  or  alternatively  P2j /c,  [30,  31]  with  Z  =  2  molecules  per  unit 
cell.  Olinger  et  al.  [27]  have  investigated  the  structure  of  the  (3-phase  within  P2(/c  symmetry 
settings  for  pressures  up  to  7.47  GPa  and  have  shown  that,  for  this  entire  pressure  range,  the 
crystal  remains  stable.  In  a  previous  study  [3],  the  thermal  expansion  properties  of  the 
crystallographic  [3-phase  described  within  P2i/n  symmetry  was  investigated.  In  order  to 
reconcile  the  two  possible  settings  of  the  same  phase  (i.e.,  with  P2[/n  or  P2i/c  symmetries),  the 
MP  and  NPT-MD  results  are  presented  for  both  of  these  symmetries.  The  corresponding  results 
are  presented  in  Table  3  and  compared  to  the  experimental  data  in  Figure  3.  The  predicted  lattice 
dimensions  and  unit  cell  volume  in  NPT-MD  simulations  at  298  K  and  zero  pressure  starting 
with  the  structure  with  P2i/n  symmetry  are  in  extremely  good  agreement  with  the  corresponding 
experimental  data.  The  percent  errors  for  lattice  dimensions  a,  b,  c,  and  unit  cell  volume  are 
-19%,  -0.72%,  0.64%,  and  0.66%,  respectively.  The  corresponding  MP  results  indicate  a 
similarly  good  agreement  with  a  maximum  deviation  of -1.79%  for  the  lattice  dimension  b  and 
-1.48%  for  the  unit  cell  volume.  Similar  NPT-MD  simulations  using  the  structure  with  the  P2]/c 
symmetry  give  results  that  are  slightly  less  accurate:  the  maximum  deviations  are  4.1%  for  lattice 
dimension  c  and  3.18%  for  the  unit  cell  volume. 

The  effect  of  increasing  the  pressure  on  lattice  dimensions  and  volume  is  also  shown  in 
Figure  3  and  Table  3.  The  results  in  the  upper  frame  correspond  to  those  obtained  from 
NPT-MD  and  MP  calculations  starting  with  the  structure  with  the  P2i/c  symmetry.  The 
predicted  dimensions  a  and  c  remain  very  close  to  the  experimental  values,  while  the  b 
dimension  deviates  slightly  with  the  increase  of  pressure.  In  Figure  3(b),  the  variation  of  the 
relative  volume  V/V0  is  compared  for  both  P2i/c  and  P2i/n  symmetries  with  corresponding 
experimental  results  [27].  The  NPT-MD  predictions  for  either  space  group  are  in  closer 
agreement  to  the  experimental  result  than  the  MP  calculations.  For  example,  the  curve  V/V0  for 
P2i/c  is  extremely  close  to  the  experimental  values  with  a  deviation  at  7.47  GPa  of  0.18%,  while, 
for  the  P2i/n  setting,  the  difference  of  1.45%.  Also,  as  seen  in  the  RDX  calculations,  MP 
predictions  for  both  lattice  dimensions  and  lattice  volume  of  (3-HMX  crystal  are  close  but  slightly 
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Table  3.  Lattice  Parameters  Obtained  in  Crystal  Packing  and  Molecular  Dynamics  Calculations  for  the  (3-HMX  Crystal  as 
a  Function  of  Pressure  at  T  =  298  K 
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Data  from  Cromer  et  al.  [29]. 

Data  from  Kohno  et  al.  [30]  and  Choi  and  Boutin  [31  ]. 
Data  from  Olinger  et  al.  [27]. 
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Figure  3.  Comparison  of  the  Crystallographic  Parameters  Obtained  in  MP  and  NPT-MD 
Simulations  for  the  fi-HMX  Crystal  With  the  Experimental  Results.  The 
Variations  of  Lattice  Dimensions  Indicated  in  the  Top  Panel  (a)  Correspond  to 
the  Structure  With  P2x/c  Symmetry  Where  Experimental  Data  Are  Available. 
The  Lower  Panel  (b)  Indicates  the  Variations  of  the  Relative  Volume  as  a 
Function  of  Pressure  for  Both  P2x/c  and  P2x/n  Symmetries  Together  With  the 
Corresponding  Experimental  Results. 


less  accurate  than  are  the  corresponding  NPT-MD  results.  This  result  is  expected  because  MP 
calculations  do  not  include  the  thermal  effects  considered  in  NPT-MD  simulations. 
Consequently,  the  predicted  lattices  are  less  compressible  than  would  be  expected  in  a 
simulation  that  correctly  incorporates  thermal  effects. 


14 


The  bulk  modulus  and  pressure  derivative  calculated  using  the  results  from  Table  3  and 
equation  (3)  are  given  in  Table  2.  The  results  are  in  reasonable  agreement  with  experimental 
values,  though  the  variation  from  experiment  is  larger  for  HMX  than  the  RDX  results.  The  P2]/c 
results  are  closer  to  experiment  than  those  determined  from  the  P2i/n  results. 

4.3  HNIW  Crystal.  HNIW,  a  polycyclic  nitramine,  has  been  characterized  as  “the  densest 
and  most  energetic  explosive  known”  [33].  It  exists  in  at  least  five  polymorphic  states,  four  of 
which  are  stable  at  ambient  conditions  and  have  been  resolved  (  a-hydrate,  e,  P,  and  y)  by  x-ray 
crystallography  [34].  The  molecular  structure  of  these  polymorphs  appears  to  be  that  of  two 
bridged  RDX  molecules  and  is  similar  to  that  shown  in  Figure  1(c)  for  the  e-phase.  The  main 
differences  between  the  configurations  of  the  different  polymorphs  are  in  the  orientation  of  the 
nitrogroups  relative  to  the  ring.  The  e  polymorph  that  is  considered  in  this  work  is  the  most 
stable  phase  at  room  temperature  [2].  It  crystallizes  in  the  P2i/n  space  group  and  has  Z  =  4 
molecules  per  unit  cell  [35]. 

The  calculated  lattice  dimensions  at  different  pressures  for  this  crystal  are  given  in  Table  4 
and  a  visual  comparison  of  the  experimental  data  of  Pinkerton  [35]  is  given  in  Figure  4.  The 
lattice  dimensions  predicted  by  MD-NPT  and  MP  simulations  are  in  very  good  agreement  with 
the  experimental  values.  For  example,  at  the  highest  pressure  considered  experimentally 
(2.5  GPa),  the  deviations  from  experiment  are  1.47%,  1.73%,  and  2.17%  for  the  a,  b,  and  c 
lattice  dimensions,  respectively.  Also,  for  this  pressure,  the  calculated  volumetric  compression 
shown  in  Figure  4(b)  is  0.90,  while  the  corresponding  experimental  value  is  0.88.  Using  the 
variation  of  the  unit  cell  volume  given  in  Table  4  with  pressure  and  equation  (3),  a  bulk  modulus 
Bo  =  15.58  GPa  and  a  pressure  derivative  Bo'  =  9.37  have  been  determined.  However,  no 
experimental  values  were  found  against  which  comparison  could  be  made  of  these  calculated 
values. 


4.4  PETN  Crystal.  The  final  system  chosen  for  assessment  of  the  interaction  potential  is 
the  non-nitramine  explosive  PETN.  This  system  could  be  considered  a  more  difficult  test  than 
the  preceding  systems,  since  the  molecular  conformation  (see  Figure  l[d])  is  much  more  floppy 
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Figure  4.  Comparison  of  the  Crystallographic  Lattice  Dimensions  Obtained  in  MP  and 
NPT-MD  Simulations  for  the  s-HNIW  Crystal  With  the  Experimental  Results. 
Dependence  of  (a)  Lattice  Dimensions  and  (b)  Volume  Compression  V/Voon  the 
External  Pressure. 

than  in  the  previous  three  systems.  This  characteristic  suggests  that  the  rigid-body 
approximation  assumed  in  the  simulations  might  be  inadequate.  Consequently,  it  is  expected 
that  predictions  of  crystallographic  parameters  for  this  type  of  crystal  within  the  constraints  of 
these  simulations  will  be  less  accurate  than  those  obtained  previously,  particularly  for  the  higher 
pressure  regime. 

The  experimental  investigations  have  shown  that  PETN  can  exist  in  two  different  phases:  a 
tetragonal  phase,  also  called  form  i  [37]  and  an  orthorhombic  phase,  known  as  form  ii  [38],  Both 
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Table  4.  Lattice  Parameters  Obtained  in  Crystal  Packing  and  Molecular  Dynamics 
Calculations  for  the  e-HNIW  Crystal  as  a  Function  of  Pressure  at  T  =  298  K 


Pressure 

(GPa) 

Molecular  Packing 

"  NPT-MD 

Lattice  Lengths 

Unit  Cell 
Volume 
(A3) 

Lattice  Lengths 

Unit  Cell 
Volume 
(A3) 

a 

o 

(A) 

b 

0 

(A) 

c 

0 

(A) 

a 

0 

(A) 

b 

0 

(A) 

c 

o 

(A) 

Exp.a 

8.8278 

12.5166 

13.3499 

1412.483 

— 

— 

— 

— 

Exp.b 

8.8847 

12.6210 

13.4310 

1439.000 

— 

— 

— 

— 

0.00 

8.8550 

12.4897 

13.4624 

1437.349 

8.8998 

12.5916 

13.5600 

1464.750 

0.50 

8.7915 

12.3701 

13.3448 

1401.346 

8.8375 

12.4457 

13.4238 

1424.500 

1.00 

8.7429 

12.2784 

13.2530 

1373.862 

8.7727 

12.3399 

13.3132 

1389.833 

1.50 

8.6696 

12.1403 

13.1119 

1351.604 

8.7440 

12.2432 

13.2245 

1365.917 

2.00 

8.6696 

12.1403 

13.1119 

1332.776 

8.6981 

12.1857 

13.1631 

1345.417 

2.50 

8.6403 

12.0851 

13.0548 

1316.499 

8.6699 

12.14465 

13.1128 

1326.120 

3.00 

8.6143 

12.0361 

13.0039 

1302.117 

8.6303 

12.06299 

13.0326 

1309.750 

3.50 

8.5911 

.11.9920 

12.9578 

1289.259 

8.6066 

12.02246 

12.9882 

1287.333 

a  Data  from  Gilardi  [35]. 


b  Data  obtained  from  Pinkerton  [36]  by  extrapolation  to  P  =  0. 


previous  experimental  results  [39]  and  theoretical  values  [5]  indicate  that  the  tetragonal  phase  is 
the  most  stable.  Therefore,  in  this  study,  the  phase  that  crystallizes  in  space  group  P421cand 
has  Z  =  2  molecules  per  unit  cell  (form  i)  is  analyzed.  The  results  of  our  MP  and  NPT-MD 
calculations  are  given  in  Table  5  and  shown  in  Figure  5. 


The  isothermal  linear  and  volume  compression  of  tetragonal  PETN  has  been  previously 
investigated  by  Olinger  et  al.  [32]  for  pressures  up  to  10  GPa  using  an  x-ray  diffraction 
technique.  The  pressure  dependence  of  the  experimental  lattice  dimensions  and  unit  cell  volume 
as  represented  in  Figure  5  were  fitted  using  a  cubic  polynomial  in  pressure  powers  given  in 
Table  2  [32].  In  the  same  table,  the  corresponding  best-fit  parameters  obtained  are  given  based 
on  predicted  NPT-MD  data  given  in  Table  5.  The  two  sets  of  fitted  parameters  have  similar 
values,  indicating  an  acceptable  agreement  between  the  experimental  and  predicted  lattice 
dimensions.  A  more  direct  comparison  is  given  in  Figure  5(a),  where  both  MD  and  MP  lattice 
dimensions  are  represented  together  with  experimental  values.  In  the  region  of  low  pressure 
(<4GPa),  the  agreement  is  very  good  with  relative  errors  of  -0.35%  and  -0.67%  for  the 
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Figure  5.  Comparison  of  the  Crystallographic  Parameters  Obtained  in  MP  and  NPT-MD 
Simulations  for  Tetragonal  Phase  of  PETN  Crystal  With  the  Experimental 
Results.  Dependence  of  (a)  Lattice  Dimensions  and  (b)  Volume  Compression 
V/Vo  on  the  External  Pressure. 

a  and  c  lattice  dimensions  and  -1.44%  for  the  unit  cell  volume.  In  addition,  as  indicated  in  Table 
5,  the  lattice  dimensions  a  and  b  remain  equal  (within  the  errors  of  simulation),  while  the  cell 
angles  (not  shown)  remain  approximately  90.0°  in  agreement  with  the  tetragonal  symmetry  of  the 
lattice.  By  increasing  the  pressure  from  0  to  9  GPa,  it  can  be  seen  in  Figure  5  that  the  predicted 
lattice  dimensions  are  very  close  to  the  experimental  values  up  to  about  6  GPa.  For  this  pressure, 
the  relative  errors  are  1.33%  and  1.25%  for  the  a  and  c  dimensions.  Above  6  GPa,  the  deviations 
of  the  predicted  values  from  the  experiment  increase  more  rapidly,  reaching  values  of  -2.08% 
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Table  5.  Lattice  Parameters  Obtained  in  Crystal  Packing  and  Molecular  Dynamics 
Calculations  for  the  PETN  Crystal  as  a  Function  of  Pressure  at  T  =  298  K 


1^ 

Molecular  Packing 

P  NPT-MD  | 

Lattice  Lengths 

Unit  Cell 

Lattice  Lengths 

Unit  Cell 

Pressure 

a 

b 

c 

Volume 

a 

b 

c 

Volume 

iraa 

\msm 

mxSM 

■II 

(A) 

0 

(A) 

(A) 

(A3) 

Exp.* 

9.3776 

9.3776 

6.7075 

589.853 

— 

— 

_ 

_ 

Exp.b 

9.3860 

9.3860 

6.7150 

591.571 

— 

— 

— 

— 

0.00 

9.2843 

9.2841 

6.5995 

568.853 

9.3467 

9.3410 

6.6546 

580.996 

0.68 

9.1950 

9.1950 

6.5214 

551.371 

9.2411 

9.2418 

6.5609 

560.414 

0.86 

9.1762 

9.1760 

6.5043 

547.667 

9.2185 

9.2202 

6.5404 

555.972 

1.28 

9.1371 

9.1370 

6.4682 

540.002 

9.1661 

9.1659 

6.4967 

545.833 

1.71 

9.1028 

9.1025 

6.4358 

533.259 

9.1275 

9.1286 

6.4603 

538.278 

2.25 

9.0654 

9.0652 

6.3997 

525.925 

9.0851 

9.0880 

6.4209 

530.139 

2.90 

9.0269 

9.0267 

6.3616 

518.363 

9.0454 

9.0450 

6.3789 

521.899 

3.42 

8.9998 

8.9997 

6.3342 

513.042 

9.0168 

9.0162 

6.3499 

516.250 

4.90 

8.9359 

8.9358 

6.2674 

500.448 

8.9476 

8.9476 

6.2791 

502.714 

5.64 

8.9093 

8.9091 

6.2384 

495.166 

8.9197 

8.9193 

6.2495 

497.199 

5.80 

8.9039 

8.9036 

6.2325 

494.092 

8.9134 

8.9132 

6.2423 

495.934 

6.73 

8.8746 

8.8743 

6.2000 

488.286 

8.8835 

8.8835 

6.2089 

489.994 

7.05 

8.8618 

8.8617 

6.1855 

485.751 

8.8722 

8.8724 

6.1983 

487.918 

7.45 

8.8540 

8.8538 

6.1767 

484.201 

8.8609 

8.8608 

6.1843 

485.559 

7.65 

8.8486 

8.8484 

6.1705 

483.125 

8.8549 

8.8554 

6.1780 

484.438 

8.11 

8.8365 

8.8364 

6.1566 

480.725 

8.8422 

8.8432 

6.1639 

481.982 

8.32 

8.8312 

8.8311 

6.1504 

479.665 

8.8376 

8.8373 

6.1571 

480.869 

8.40 

— 

— 

— 

— 

8.8355 

8.8354 

6.1543 

480.441 

9.04 

— 

— 

— 

— 

8.8195 

8.8192 

6.1357 

477.250 

a  Data  from  Trotter  [37]. 
b  Data  from  Olinger  et  al.  [32]. 


and  1 .26%  for  the  a  and  c  dimensions,  respectively,  at  P  =  8.84  GPa.  This  trend  is  accentuated  in 
the  variation  of  the  relative  volume  V/V0,  where  the  deviation  at  5  GPa  is  about  4.6%,  while,  at 
8.84  GPa,  the  deviation  reaches  a  value  of  7.8%.  The  lower  compressibility  of  the  theoretical 
model  is  also  reflected  in  the  bulk  modulus.  The  predicted  bulk  modulus  using  equation  (3)  and 
the  data  in  Table  5  is  14.0  GPa,  while  the  experimental  value  is  9.9  GPa  [32].  Thus,  it  is  evident 
that  the  assumption  of  rigid  molecules  in  these  molecular  simulations  is  invalid  for  pressures 
above  5  GPa. 
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Changes  in  volume  with  pressure  upon  uniaxial  compression  in  either  dimension  have  also 
been  explored  using  this  model.  This  investigation  is  accomplished  using  the  pressure 
dependencies  of  the  individual  lattice  parameters  determined  in  Table  2,  which  are  assumed  to 
remain  valid  in  the  case  of  uniaxial  compression.  The  results  presented  in  Figure  6  clearly 
indicate  that  for  a  given  applied  external  pressure  the  variation  of  the  volume  is  larger  in  the  case 
of  compression  along  c-axis  than  in  the  case  of  a-axis.  This  result  indicates  a  strong  sensitivity 
of  the  volume  change  with  the  type  of  compression  applied  and  is  in  agreement  with  the  findings 
reported  by  Kunz  [40]  based  on  ab  initio  periodic  Hartree-Fock  calculations. 


Volume  (A3) 

Figure  6.  The  Pressure- Volume  Dependence  for  PETN  in  the  Case  of  Uniaxial 
Compression  Along  a-Axis  and  c-Axis,  Respectively.  The  Calculated  Curves 
Have  Been  Determined  Using  the  Pressure  Coefficients  Indicated  in  Table  2. 


Considering  the  ensemble  of  results  presented  in  this  work,  it  appears  that  the  rigid-molecule 
model  can  be  used  in  connection  with  the  present  intermolecular  potential  within  a  carefully 
considered  range  of  pressures.  Depending  on  the  type  of  molecular  structure  (i.e.,  more  floppy 
as  PETN,  or  more  rigid  as  [poly]cyclic  nitramines  RDX,  HMX,  and  HNIW),  the  range  of 
pressures  in  which  realistic  predictions  can  be  made  varies  from  about  5  GPa  in  the  case  of 
PETN  to  more  than  7.5  GPa  for  HMX.  The  main  effect  of  the  initial  compression  of  these 
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crystals  is  represented  by  reduction  of  the  intermolecular  distances.  Such  an  effect  can  be  seen 
in  Figure  7,  where  the  RDFs  of  the  COM-COM  pairs  for  all  molecules  in  the  MD  simulation  box 
are  represented.  The  main  changes  in  the  RDF  distributions  are  the  shifts  of  the  peaks  toward 
smaller  distance  values  with  increasing  pressure,  indicating  the  compression  of  these  materials. 
In  addition,  there  is  a  continuous  decrease  of  the  amplitude  of  molecular  oscillations  around  the 
equilibrium  positions  reflected  by  additional  resolved  peaks,  particularly  for  the  HMX  and 
HNIW  structures. 


5.  Conclusions 

The  hydrostatic  compression  effects  have  been  investigated  on  some  important  energetic 
materials  RDX,  HMX,  HNIW,  and  PETN  through  crystal  packing  and  isothermal-isobaric 
molecular  dynamics  simulations.  The  potential  used  in  these  calculations  was  previously 
developed  for  RDX  and  shown  to  be  transferable  to  30  nitramine  crystals  [4]  and  to  51  other 
non-nitramine  crystals  [5]  consisting  of  molecules  that  contained  functional  groups  associated 
with  energetic  materials.  These  systems  include  different  types  of  nitroalkanes,  nitroaromatics, 
nitrocubanes,  polynitroadamantes,  hydroxy-nitroderivatives,  nitrobenzonitriles,  nitrobenzotri/ 
azoles,  and  nitrate  esters. 

The  tests  of  this  potential  indicate  that  the  predictions  of  the  crystallographic  parameters  for 
RDX,  HMX,  and  HNIW  are  very  good  with  lattice  dimension  errors  below  2.2%  for  RDX,  4.1% 
for  HMX,  and  2.17%  for  HNIW.  Moreover,  this  potential  is  able  to  predict  both  the  changes  of 
the  crystallographic  structures  with  pressure  and  the  bulk  moduli  and  pressure  derivatives.  For 
RDX,  the  predicted  bulk  modulus  at  zero  pressure  is  Bo  =  12.93  GPa,  while  the  experimental 
value  is  BoeXp  =13.0  GPa.  Similarly,  for  HMX  Bo  =  14.64  GPa  vs.  an  experimental  value  Boexp  = 
13.5  GPa. 
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Figure  7.  Radial  Distribution  Functions  for  COM-COM  Pairs  as  Function  of  Pressure  for 
(a)  a-RDX,  (b)  p-HMX,  (c)  e-HNIW;  and(d)  PETN  (Tetragonal  Phase). 

Another  finding  of  the  present  work  is  that  molecular-packing  calculations  give 
crystallographic  parameters  that  are  very  close  to  those  determined  in  NPT-MD  simulations  but 
at  a  fraction  of  cost  in  the  computing  time.  This  result  indicates  the  utility  of  MP  calculations 
not  only  for  equilibrium  crystallographic  structures  but  also  when  the  compression  effects  are  of 
interest. 

In  the  case  of  crystals  with  floppy  molecules  such  as  PETN,  the  range  of  pressure  over  which 
the  present  model  can  make  accurate  predictions  is  more  limited,  in  this  case,  to  about  5  GPa. 
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This  result  does  not  represent  a  limitation  of  the  potential  parameters  involved  but  rather 
indicates  that  the  rigid  molecular  model  used  fails  to  be  valid  for  such  systems. 

The  success  of  the  present  potential  energy  parameters  in  describing  not  only  the  ambient 
equilibrium  structures  of  different  crystals  with  functional  groups  associated  with  explosives  but 
also  the  effects  of  external  stimuli,  such  as  pressure  and  temperature,  provides  significant 
incentive  to  further  develop  this  model  by  incorporating  the  intramolecular  degree  of  freedom. 
This  will  be  done  in  future  work. 
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